Automatic finite element implementation of hyperelastic material with a double numerical differentiation algorithm

In order to accelerate implementation of hyperelastic materials for finite element analysis, we developed an automatic numerical algorithm that only requires the strain energy function. This saves the effort on analytical derivation and coding of stress and tangent modulus, which is time-consuming and prone to human errors. Using the one-sided Newton difference quotients, the proposed algorithm first perturbs deformation gradients and calculate the difference on strain energy to approximate stress. Then, we perturb again to get difference in stress to approximate tangent modulus. Accuracy of the approximations were evaluated across the perturbation parameter space, where we find the optimal amount of perturbation being 10 to obtain stress and 10 to obtain tangent modulus. Single element verification in ABAQUS with Neo-Hookean material resulted in a small stress error of only 7 × 10 on average across uniaxial compression and tension, biaxial tension and simple shear situations. A full 3D model with Holzapfel anisotropic material for artery inflation generated a small relative error of 4 × 10 for 1 Corresponding author. 2 inflated radius at 25 kPa pressure. Results of the verification tests suggest that the proposed numerical method has good accuracy and convergence performance, therefore a good material implementation algorithm in small scale models and a useful debugging tool for large scale models.